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Coupling local, slowly adapting variables to an attractor network allows to destabilize all attractors, turning 
them into attractor ruins. The resulting attractor relict network may show ongoing autonomous latching 
dynamics. We propose to use two generating functionals for the construction of attractor relict networks, a 
Hopfield energy functional generating a neural attractor network and a functional based on 
information-theoretical principles, encoding the information content of the neural firing statistics, which 
induces latching transition from one transiently stable attractor ruin to the next. We investigate the 
influence of stress, in terms of conflicting optimization targets, on the resulting dynamics. Objective 
function stress is absent when the target level for the mean of neural activities is identical for the two 
generating functionals and the resulting latching dynamics is then found to be regular. Objective function 
stress is present when the respective target activity levels differ, inducing intermittent bursting latching 
dynamics. 

The use of objective functions for the formulation of complex systems has seen a study surge of interest. 
Objective functions, in particular objective functions based on information theoretical principles 1 " 7 , are 
used increasingly as generating functionals for the construction of complex dynamical and cognitive 
systems. There is then no need to formulate by hand equations of motion, just as it is possible, in analogy, 
to generate in classical mechanics Newton's equation of motion from an appropriate Lagrange function. When 
studying dynamical systems generated from objective functions encoding general principles, one may expect to 
obtain a deeper understanding of the resulting behavior. The kind of generating functional employed also serves, 
in addition, to characterize the class of dynamical systems for which the results obtained may be generically 
valid. 

Here we study the interplay between two generating functionals. The first generating functional is a simple 
energy functional. Minimizing this objective function one generates a neural network with predefined point 
attractors, the Hopfield net 8 ' 9 . The second generating functional describes the information content of the indi- 
vidual neural firing rates. Minimizing this functional results in maximizing the information entropy 6 ' 7 and in the 
generation of adaption rules for the intrinsic neural parameters, the threshold and the gain. This principle has 
been denoted polyhomeostatic optimization 710 , as it involves the optimization of an entire function, the distri- 
bution function of the time-averaged neural activities. 

We show that polyhomeostatic optimization destabilizes all attractors of the Hopfield net, turning them into 
attractor ruins. The resulting dynamical network is an attractor relict network and the dynamics involves 
sequences of continuously latching transient states. This dynamical state is characterized by trajectories slowing 
down close to the succession of attractor ruins visited consecutively. The two generating functionals can have 
incompatible objectives. Each generating functional, on its own, would lead to dynamical states with certain 
average levels of activity. Stress is induced when these two target mean activity levels differ. We find that the 
system responds to objective function stress by resorting to intermittent bursting, with laminar flow interseeded 
by burst of transient state latching. 

Dynamical systems functionally equivalent to attractor relict networks have been used widely to formulate 
dynamics involving on-going sequences of transient states. Latching dynamics has been studied in the context of 
the grammar generation with infinite recursion 1112 and in the context of reliable sequence generation 1314 . 
Transient state latching has also been observed in the brain 1516 and may constitute an important component 
of the internal brain dynamics. This internal brain dynamics is autonomous, ongoing being modulated, but not 
driven, by the sensory input 1718 . In this context an attractor relict network has been used to model autonomous 
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neural dynamics in terms of sequences of alternating neural firing 
patterns 19 . The modulation of this type of internal latching dynamics 
by a stream of sensory inputs results, via unsupervised learning, in a 
mapping of objects present in the sensory input stream to the pre- 
existing attractor ruins of the attractor relict network 20,21 , the asso- 
ciative latching dynamics thus acquiring semantic content. 

Results 

As an introductory example we consider one of the simplest possible 
networks having two distinct attractors in terms of minima of the 
energy functional, as defined by Eq. (7). The dynamics of the 3-site 
network illustrated in Fig. 1 is given by 



k\ = — Txi + w + y 2 — w y 3 
x 2 = -Tx 2 + w + (y l +y 3 ) , 
x 3 = — Tx 3 + w + y 2 — w~yi 



(i) 



with w + > 0 and w~ > 0 denoting the excitatory and inhibitory link 
strength respectively. Here the xjy { correspond to neural membrane 
potentials and firing rates, respectively, related through a sigmoidal 
transfer function, see Eq. (5), which is parameterized by the gain a { 
and the threshold b { . For fixed intrinsic parameters a = a t and b = b b 
viz for vanishing adaption rates e fl , = 0, the network has two pos- 
sible phases, as shown in Fig. 1. There is either a single global attrac- 
tor, with activities y { determined mostly by the value of the threshold 
b, or two stable attractors, with either y 1 or y 3 being small and the 
other two firing rates large. 

The point attractors present for vanishing intrinsic adaption 
become unstable for finite adaption rates e fl , >0. Point attractors, 
defined by Jc f - = 0 = y-, have time-independent neural activities. The 
objective (9), to attain a minimal Kullback-Leibler divergence, can 
however not be achieved for constant firing rates y { . The polyhomeo- 
static adaption (11) hence forces the system to become autono- 
mously active, as seen in Fig. 2. 

In Fig. 2 we present the results for the overlaps O p and A p , as 
defined by Eqs. (13) and (14), of the neural activity with the stored 
memories. One observes, that the original attractors become unstable 
in the presence of finite polyhomeostatic adaption, retaining how- 
ever a prominent role in phase space. Unstable attractors, transiently 
attracting the phase space flow, can be considered to act as 'attractor 
relicts'. The resulting dynamical system is hence an attractor relict 



network, viz a network of coupled attractor ruins 19 ' 20 . When coupling 
an attractor network to slow intrinsic variables, here the threshold b 
and the gain a, the attractors are destroyed but retain presence in the 
flow in terms of attractor ruins, with the dynamics slowing down 
close to the attractor relict. 

The concept of an attractor ruin is quite general and attractor relict 
networks have implicitly been generated in the past using a range of 
functionally equivalent schemes. It is possible to introduce additional 
local slow variables, like a reservoir, which are slowly depleted when a 
unit is active 19 ' 20 . Similarly, slowly accumulating spin variables have 
been coupled to local dynamic thresholds in order to destabilize 
attractors 13 , and local, slowly adapting, zero state fields have been 
considered in the context of latching Potts attractor networks 1112 ' 22 ' 23 . 
Alternatively, attractors have been destabilized in neural networks by 
introducing slowly adapting asymmetric components to the inter- 
neural synaptic weights 14 . 

Phase boundary adaption. In Fig. 3 we have superimposed, for the 
three- site network, a typical set of trajectories, with finite adaption 
rates e a = 0.1 and q, = 0.01, onto the phase diagram evaluated for 
non-adapting neurons, with e a = 0 = as shown in Fig. 1. The three 
trajectories (green: (a 2 , b 2 ), red and blue: (a u b x ) and (a 3 , b 3 )) all start 
well in the region with a single global fixed point, at ^ = 5, b t = — 0.5 
(i = 1, 2, 3). Two features of the flow are eye-catching. 

• The intrinsic parameters of the active neurons settle, after a tran- 
sient initial period, at the phase boundary and oscillate alternating 
between the phase with a single stable fixed point and the phase 
with two attractors. 

• The trajectory (a 2 , b 2 ) for the central neuron adapts to a large 
value for the threshold, taking care not to cross any phase bound- 
ary during the initial transient trajectory. 

We first discuss the physics behind the second phenomenon. The 
Euclidean distance in phase space between the two attractors, when 
existing, can be used as an order parameter. In Fig. 4 the order 
parameter is given for vertical cuts through the phase diagram. 
The lower transition is evidently of second order, which can be con- 
firmed by a stability analysis of (1), with the upper transition being of 
first order, in agreement with a graphical analysis of (1). A first order 
transition will not be crossed, in general, by an adaptive process, as 
there are no gradients for the adaption to follow. This is the reason 
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Figure 1 | Left: A three-site graph with symmetric excitatory (solid green lines) and inhibitory connections (red dashed lines). Right: The corresponding 
phase diagram, for fixed T = 1 and w + = 1 = w~, as function of the gain a = ai and threshold b = b iy for i = 1, 2, 3. For large/small thresholds the sites tend 
to be inactive/active. The color encodes the activity of the most active neuron. Inside the shaded area there are two stable attractors, outside a single, 
globally attracting state. The binary representations of the two non-trivial attractors are (1,1,0) and (0,1,1), where 1/0 denotes an active/inactive site. The 
lower and upper lines (blue/white) of the phase transition are of second and first order respectively, with the lower phase transition line determined by 
1 = ay{\ —y) where y=y\ =y 3 is the activity of the single fixed point at site one and three. The two critical lines meet at (a a b c ) = (4, 0.413) (blue dot). 
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Figure 2 | Time series of the three-site network shown in Fig. 1, for T = l,e a = 0.1, ey = 0.01, Ai = 0>A 2 = 0 and wr = 1. The dynamics is given by (6) 
and (11). The dynamics retraces periodically the original attractor states, as one can see from the oscillation of the overlap O p and the A py compare Eqs. 
(13) and (14), between the patterns of neural and attractor activities. 



that the trajectory for (a 2 , b 2 ) shown in Fig. 3 needs to take a large 
detour in order to arrive to its target spot in phase space. This beha- 
vior is independent of the choice of initial conditions. 

Using an elementary stability analysis one can show that the con- 
dition for the second-order line is 1 = ay(l —y), where y=y\ =y^ 



(for the single global fixed point). Only a single global attractor is 
consequently stable for a < 4 (since ye[0 9 1] and j>(l —y) <0.25) and 
the second and the first order lines meet at (a c , b c ) = (4, 0.413), where 
the critical threshold b c is determined by the self- consistent solution 
of b c = 1/[1 + exp(4b c — 4)] — 1/2. The trajectories (a l5 b x ) and 




Figure 3 | The trajectories of (fl;(t), bi(t)), for finite adaption rates, e a = 0.1 and = 0.01, superimposed onto the phase diagram of the three-site 
network shown in Fig. 1. The color encodes the activity of the most active site. Inside the shaded blue area two attractors are stable, outside there is only a 
single global attractor. The parameters are the same as for the simulation presented in Fig. 2. The thresholds bi(t) and b 3 (t) (red/blue trajectories) oscillate 
across the second-order phase boundary (blue line), the threshold b 2 (t) (green trajectory) acquires a large value, receiving two excitatory inputs, and 
avoids the first-order phase boundary (white line) during the initial transient. All three trajectories start in the center, lower-half of the phase diagram, at a 
= 5, b = —0.5 (black dot). Note that the underlying phase diagram is for identical thresholds bj and gains a z . 
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Figure 4 | Vertical cuts (for constant gains a = 4, 0, 4.3, 4.7, 5.0, 6.0, 
bottom to top curves) through the phase diagram of the 3-site graph shown 



in Fig. 1. Plotted is the order parameter y^ j=1 yy\ 



,(!) 



-yf^ , where 



are the two fixed point solutions (y = 1, 2) in the 
coexistent region. The lower/upper transitions are second/first order. 

(a 3 , b 3 ) are hence oscillating across the locus in phase space where 
there would be a second-order phase transition, compare Fig. 3, for 
the case of identical internal parameters a t and b { . With adaption, 
however, the respective internal parameters take distinct values, (a, b 
± db) for the first/third neuron and (a, b 2 ) for the second neuron, 
with a ~ 6, b ~ 0 and fr 2 ~ 1- In fact the system adapts the intrinsic 
parameters dynamical in a way that a non- stopping sequence of 
states 



(2) 



is visited, where f 1 = (1, 1, 0), = (0, 1, 1) denote the two stable 
cliques and £° = (1, 1, 1) the fixed point in the phase having only a 
single global attractor. It is hence not a coincidence that the system 
adapts autonomously, as shown in Fig. 3, to a region in phase space 
where there would be a second order phase transition for identical 
internal parameters. At this point, small adaptive variations make the 
sequence (2) achievable. The adaption process hence shares phe- 
nomenological some aspects with what one calls 'self-organized 
quasi-criticality' in the context of absorbing phase transitions 27 ' 28 , 
which denotes the circumstance that a system without energy con- 
servation may become critical only when tuning an external rechar- 
ging rate. For imperfect tuning the trajectories would linger close to 
second-order phase transition point, without actually reaching it, a 
phenomenology close to one seen in Fig. 3. 

Objective function stress. We will now investigate networks of larger 
size N. In principle we could generate random synaptic link matrices 
{Wfj} and find their eigenstates f = ({?,... ,{^) using a diagonali- 
zation routine. When using the Hopfield encoding (15) for the 
synaptic weights Wy the attractors are known, corresponding, as 



long as the number N p of patterns is not too large 8 ' 9 , to the to the 
stored patterns We here make use of the Hopfield encoding for 
convenience, no claim is made that memories in the brain are actually 
stored and defined by (15). 

We consider in the following random binary patterns, as illu- 
strated in Fig. 5, 

^ p ( 1 with probability a 
1 \ 0 with probability 1 — a 

where a is the mean activity level or sparseness. The patterns have 
in general a finite, albeit small, overlap, as illustrated in Fig. 5. The 
target distribution q(y) for the intrinsic adaption has, for k 2 — 0, an 
expected mean 



fdy 
Jo 



yq(y) = \- 



(4) 



which can be evaluated noting that the support of q(y) is [0, 1]. The 
target mean activity fi can now differ from the average activity of an 
attractor relict, the mean pattern activity a, compare (3). The differ- 
ence between the two quantities, viz between the two objectives, 
energy minimization vs. polyhomeostatic optimization, induces 
stress into the latching dynamics. 

In Fig. 6 we present the time evolution of the overlaps A p and O p 
for a = 0.3 = fi and adaption rates e a = 0.1, e\, = 0.01. In this case the 
target mean activity, fi, of the intrinsic adaption rules (1 1) is consist- 
ent with the mean activity a of the stored patterns £ p , viz with the 
mean activity level of the attractor relicts. One observes that the 
system settles into a limiting cycle with all seven stored patterns 
becoming successively transiently active, a near- to -perfect latching 
dynamics. The dynamics is very stable and independent of initial 
conditions, which were selected randomly. 

In Fig. 7 we present the time evolution, for 20 out of the N = 100 
neurons, of the respective individual variables. The simulation para- 
meters are identical for Figs. 6 and 7. Shown in Fig. 7 are the indi- 
vidual membrane potentials #,-(f)> the firing rates y^t), the gains a { {t) 
and the thresholds bj(t). The latching activation of the attractor 
relicts seen in Fig. 6 reflects in corresponding transient activations 
of the respective membrane potentials and firing rates. The oscilla- 
tions in the thresholds b^t) drive the latching dynamics, interest- 
ingly, even though the adaption rate is larger for the gain. 

The synaptic weights are symmetric and consequently also the 
overlap matrix presented in Fig. 5. The latching transitions evident in 
Fig. 6 are hence spontaneous in the sense that they are not induced by 
asymmetries in the weight matrix 14 . We have selected uncorrelated 
patterns and the chronological order of the transient states is hence 
determined by small stochastic differences in the pattern overlaps. It 
would however be possible to consider correlated activity patters £ p 
incorporating a rudimental grammatical structure 22,23 , which is how- 
ever beyond the scope of the present study. 

Stress-induced intermittency. In Fig. 8 we present the time evolu- 
tion of the overlaps A p and O p for the case when the two objective 
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Figure 5 | Left: Example of N P = 7 random binary patterns compare (3), with sparseness a = 0.3 for a network of N = 100 neurons. Right: Respective 
pairwise mutual pattern overlaps 0(£ p , £ q ). Patterns maximally overlap with themselves, the inter-pattern overlap is random and small. 



SCIENTIFIC REPORTS | 3 : 2042 | DOI: 1 0.1 038/srep02042 



4 




1000 



1100 



1200 



1300 



1400 



1500 




1000 



1100 



1200 



1300 



t 



1400 



1500 



Figure 6 | Overlaps O p and A p (color encoded, see Eqs. (13) and ( 14)), for a N = 100 network, of the neural activities with the N p = 7 attractor ruins, as a 
function of time t. The sparseness of the binary patterns defining the attractor ruins is a = 0.3. Here the average activity a of the attractor relicts and the 
target mean activity level fi = 0.3 have been selected to be equal, resulting in a limiting cycle with clean latching dynamics. 



functions, the energy functional and the polyhomeostatic optimiza- 
tion, incorporate conflicting targets. We retain the average spar- 
seness a = 0.3 for the stored patterns, as for Fig. 6, but reduced the 
target mean firing rate to fi = 0.15. This discrepancy between a and fi 
induces stress into the dynamics. The pure latching dynamics, as 
previously observed in Fig. 6, corresponds to a mean activity of 
about 0.3, in conflict with the target value fi = 0.15. 

Phases of laminar flow are induced by the objective function stress, 
and the latching dynamics occurs now in the form of intermittent 
bursts. The neural activity, see A p in Fig. 8, is substantially reduced 



during the laminar flow and the time averaged mean firing rate such 
reduced towards the target activity level of fi = 0.15. The trajectory 
does not come close to any particular attractor relict during the 
laminar flow, the overall O p being close to 0.5 for all N p = 7 stored 
pattern. 

The time evolution is hence segmented into two distinct phases, a 
slow laminar phase far from any attractor and intermittent phases of 
bursting activity in which the trajectories linger transiently close to 
attractor relicts with comparatively fast (relative to the time scale 
of the laminar phase) latching transitions. Intermittent bursting 
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Figure 7 | Time evolution, for a selection of 20 out of AT = 100 neurons, of the membrane potential x iy gain a iy threshold bj and firing rate y z - for the time 
series of overlaps presented in Fig. 6. 
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Figure 8 | Time evolution of overlaps O p and A p (color coded, compare Eqs. (13) and (14)) for all N P = 7 binary patterns, with sparseness a = 0.3, 
for a network of N = 100 neurons. The target mean neural activity is /u = 0.15, the difference between the two objective functions, viz between /u and a 
induces stress. The intermittent latching dynamics has a mean activity of about 0.3, which is too large. The phases of laminar flows between the 
burst of latching has a reduced average activity, thus reducing the time-averaged mean activity level toward the target fi = 0.15. 



dynamics has been observed previously in polyhomeostatically adap- 
ting neural networks with randomly selected synaptic weights 710 , the 
underlying causes had however not been clear. The results presented 
above show that the concept of competing generating functionals can 
be used to define objective function stress and that a given difference 
in target mean activities is the underlying force driving the system 
into the intermittent bursting regime. 

Robustness of transient state dynamics. The latching dynamics 
presented in Figs. 6 and 8 is robust with respect to system size N. 
We did run the simulation for different sizes of networks with up to N 
= 10 5 neurons, a series of sparseness parameters a and number of 
stored patterns N P . As an example we present in Fig. 9 the overlap O p 
for N = 1000 neurons and N P = 20 binary patterns with a sparseness 
of a = 0.2. No stress is present, the target mean activity level is [i = 
0.2. The latching dynamics is regular, no intermittent bursting is 
observed. There is no constraint, generically, to force the limiting 
cycle to incorporate all attractor relicts. Indeed, for the transient state 
dynamics presented in Fig. 9, some of the stored patterns are never 
activated. 

The data presented in Figs. 6, 8 and 9 is for small numbers of 
attractor relicts N p , relative to the systems size N, a systematic study 
for large values of N p is beyond the scope of the present study. 
Latching dynamics tends to break down, generically speaking, when 
the overlap between distinct attractors, as shown in Fig. 5, becomes 
large. The autonomous dynamics then becomes irregular. 

Discussion 

The use of generation functionals has a long tradition in physics in 
general and in classical mechanics in particular. Here we point out 
that using several multivariate objective functions may lead to novel 
dynamical behaviors and an improved understanding of complex 
systems in general. We propose in particular to employ generating 
functionals which are multivariate in the sense that they are used to 
derive the equations of motion for distinct, non- overlapping subsets 
of dynamical variables. In the present work we have studied a neural 
network with fast primary variables #,-(f)> the membrane potentials, 
and slow secondary variables a f -(f) and &,■(£), characterizing the 
internal behavior of individual neurons, here the gains a* and the 



thresholds b { . The time evolution of these sets of interdependent 
variables is determined respectively by two generating functionals. 

The equations of motion (leaky integrator) for the primary 
dynamical variables kj(t), the individual membrane potentials, are 
generated minimizing an energy functional. The Kullback-Leibler 
divergence is, on the other side, an example of an information- 
theoretical fucntional. Minimizing the Kullback-Leibler divergence 
between the distribution p^y) of the time-average neural firing rate)// 
and a target distribution function q(y) maximizing information 
entropy generates intrinsic adaption rules a and b for the gain a 
and the threshold b (polyhomeostatic optimization). 

Generating functionals may incorporate certain targets or con- 
straints, either explicitly or implicitly. We denote the interplay 
between distinct objectives incorporated by competing generating 
functionals "objective functions stress". For the two generating func- 
tionals considered in this study there are two types of objective func- 
tions stress. 

The minima of the energy functional are time-independent point 
attractors leading to firing-rate distributions pi(y) which are sharply 
peaked. The target firing-rated distribution q(y) for the information- 
theoretical functional is however smooth (polyhomeostasis). This 
functional stress leads to the formation of an attractor relict network. 
The mean target neural firing rate fi is, on the other side, a scalar 
parameter for the target firing-rate distribution q(y), and hence 
encoded explicitly within the information theoretical functional. 
The local minima of the energy functional, determined by the syn- 
aptic weights Wfj, determine implicitly the mean activity levels a of 
the corresponding point attractors. Scalar objective function stress is 
present for a ^ /n. 

For the two generating functionals considered, we find that the 
scalar objective function stress induces a novel dynamical state, char- 
acterized by periods of slow laminar flow interseeded by bursts of 
rapid latching transitions. We propose that objective function stress 
is a powerful tool, in general, for controlling the behavior of complex 
dynamical systems. The interplay between distinct objective func- 
tions may hence serve as a mechanism for guiding self organiza- 
tion 30 " 32 , a process also denoted targeted self organization 25 . 

The main focus of this study has been to to investigate and to quan- 
tify, for larger size networks, the competition between an objective 
function based on energy minimization and a generating functional 
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Figure 9 | Time evolution, for N = 1000, of the overlaps O p for all N P = 20 binary patterns (vertically displaced) with sparseness a = 0.2 and a target mean 
activity of fi = 0.2. There is no objective-function stress and the latching is regular, the adaption rates are e a = 0.1, = 0.01. 



based on the principle of homeostasis. Within this investigation we 
also detailed out the phase diagram of a model 3 -site network having 
either one or two transiently stable attractors. This application to a 
simple model network allowed to investigate the adaptive flow within 
the underlying phase diagram, as generated in the adiabatic limit by 
the energy functional, compare Fig. 3. The results indicate that poly- 
homeostatic adaption has the tendency to drive the system toward 
regions of phase space close to a second order phase transition, as this 
offers to generate the desired variability in neural activity. We plan to 
investigate further this interesting phenomenon in future studies. 

Methods 

We use N rate encoding neurons in continuous time, with a non-linear transfer 
function, 



where Xi(t)eTZ is the membrane potential and y ; (f) e [0, 1] the firing rate of the zth 
neuron. The dynamics of the neural activity is 

N 

Xi=-Txi+^Wijy h (6) 
;=i 

which describes a leaky integrator 24 , with Y being the leak rate and the the inter- 
neural synaptic weights. This kind of dynamics can be derived minimizing a simple 
energy functional 

i ij 

with respect to the membrane potential x, 8,9 . When using an energy functional to 
derive (6), the resulting synaptic weights are necessarily symmetrized, viz Wy = Wp. 
Here we do indeed consider only symmetric synaptic links, the dynamics of a network 
of leaky integrators would however remain well behaved when generalizing to non- 
symmetric weight matrices. 

The transfer function (5) contains two intrinsic parameters, a and b, which can be 
either set by hand or autonomously adapted, an approach considered here. The 
intrinsic parameters determine the shape of the individual firing-rated distributions 
pi(y), given generically by 

p i (y)= l -^5(y-y i (t-t))dx, (8) 

where S(') is the Dirac (5-function, and T—> °o the observation time-interval. The 
Kullback-Leibler divergence 25 



D KL (j, h q)=^dy Pi (y)lo g ^ (9) 

measures the distance between the neural firing rate distribution p t (y) and a nor- 
malized distribution q(y), with D KL ^ 0 generically and D KL = 0 only for pi(y) = q(y). 
Gaussian distributions 

q{y)ccexp(A,i y + A 2 f) (10) 

maximize the information entropy whenever mean and standard deviation are 
given 25 . Minimizing the Kullback-Leibler divergence (9) with respect to the individual 
intrinsic parameters a f and b { is hence equivalent to optimizing the information 
content of the neural activity 6 . The optimization can be performed using variational 
calculus 7,10 ' 26 , one obtains 

ai = e a (l/ai + (xi-bi)6) 

(11) 

with 6 = 1 — 2y { + (X x + 2X 2 y\) ( 1 — yd Ji and where the e a and Cb are adaption rates. 
When these adaption rates are small, with respect to the time scale of the cognitive 
dynamics (6), an implicit time average with respect to the distribution of input signals 
is performed 10,26 , compare (8), and the respective firing rate distribution p t {y) 
approaches the target distribution function q(y). The rules (11) are also denoted 
stochastic adaption rules, since they correspond 7,10,26 , for small adaption rates e a and 
££, to an implicit (stochastic) average over all input signals x t . The adaption rules (11) 
implement the optimization of an entire distribution function, and not of a single 
scalar quantity, and are hence equivalent to a polyhomeostatic optimization process 7 . 

The original network (6) has point attractors, as given by jc,- = 0, Vz. These attractors 
are destabilized for any e a ,Cb> 0, the coupled system (6) and (11) has no stationary 
points, as defined by x = 0 = a = b. Stationary point attractors lead to neural firing 
statistics far from the target distribution (10) and to a large Kullback-Leibler diver- 
gence D KL , and hence to finite stochastic gradients (11). The principle of polyho- 
meostatic adaption is hence intrinsically destabilizing. 

In deriving the stochastic adaption rules (11) one rewrites the Kullback-Leibler 
divergence D KL as an integral over the distribution of the respective membrane 
potential x b namely as D KL = J dx pi(x)d(x), with an appropriate kernel d(x). The 
optimal overall adaption rules depend on the specific shape of pi(x). The optimal 
adaption rates not dependent on the distribution of the membrane potential are, on 
the other side, given by minimizing the kernel, accdd(x)/da and bccdd(x)/db 
respectively, which leads to the adaption rates (11), which are instantaneous in 
time 7,10,26 . 

For the analysis of the simulation results we use with 



?=(*?,... ,4), wa-Jt^) 2 (12) 

the binary activity patterns of the attractors of (6). The are ^ = (1, 1, 0) and £ 2 = (0, 1, 
1 ) for the case of the N = 3 network presented in Fig. 1 . We use two criteria in order to 



Table 1 | Relation between the parameter Ai, for A 2 = 0, and the mean value \i, as given by Eq. 4, for the target distribution g(y) 
\i 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

h -9.995 -4.801 -2.672 -1.229 0 1.229 2.672 4.801 9.995 
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investigate to which degree the original binary attractors are retraced. The first cri- 
terion is the overlap O p <e [0,1], 

ii p ,y) 



°P = Y^> <CV)=E^. (13) 



which corresponds to the cosine of the angle between the actual neural activity (y x , 
y N ) and the attractor state As a second measure, of how close the actual activity 
pattern and the original attractors are, we consider the reweighted scalar product A p 
e [0, 1] 

Note that the are positive, representing the neural firing rate in attractor states. For 
the general network analysis we use the Hopfield encoding 8,9 

H*«^^£(tf-^te-5) (is) 



'a(N-l) 



P =i 



for the synaptic weights w t j, where the ^ = — $ are the arithmetic means of 

all pattern activities, for the respective sites. Here N p is the number of encoded 
attractors and a e [0, 1] the mean pattern activity, 



N 



(16) 



When using the Hopfield encoding (15), the attractors are known to correspond, to a 
close degree, to the stored patterns as long as the number N p of patterns is not too 
large 8 ' 9 . 

We simulated Eqs. (6) and (11) using 4th order classical Runge-Kutta 29 and a 
timestep of At = 0.1. The resulting dynamics is dependent on the magnitude of the 
adaption rates for the gain a and for the threshold b. In general latching dynamics is 
favored for small e\, ~ 0.01 and somewhat larger e a ~ 0.01 ... 1. We kept a constant 
leak rate T = 1. The results presented in Figs. 6 and 7 have been obtained setting X 2 = 
0 and using Table 1 for X\. 
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